Analytical probabilistic approach to the ground 
state of lattice quantum systems: exact results in 
terms of a cumulant expansion 

Massimo Ostilli^^ and Carlo Presilla^^ '^ 

^ Dipartimento di Fisica, Univcrsita di Roma "La Sapionza" , Piazzalo A. Moro 
2, Roma 00185, Italy 

^ Center for Statistical Mechanics and Complexity, Istituto Nazionalc per la 
Fisica della Materia, Unita di Roma 1, Roma 00185, Italy 

^ Istituto Nazionalc di Fisica Nucleare, Sezione di Roma 1, Roma 00185, Italy 

Abstract. We present a large deviation analysis of a recently proposed 
probabilistic approach to the study of the ground-state properties of lattice 
quantum systems. The ground-state energy, as well as the correlation functions 
in the ground state, are exactly determined as a series expansion in the cumulants 
of the multiplicities of the potential and hopping energies assumed by the system 
during its long-time evolution. Once these cumulants are known, even at a finite 
order, our approach provides the ground state analytically as a function of the 
Hamiltonian parameters. A scenario of possible applications of this analyticity 
property is discussed. 
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1. Introduction 

The real- or imaginary-time dynamics of systems described by a finite Hamiltonian 
matrix, representing bosonic or fermionic degrees of freedom, admits an exact 
probabilistic representation in terms of a proper collection of independent Poisson 
processes For a lattice system, the Poisson processes are associated to the 

links of the lattice and the probabilistic representation leads to an optimal algorithm 
which coincides with the Green Function Quantum Monte Carlo method in the limit 
when the latter becomes exact 0]. 

Recently, in we have exploited the above probabilistic representation to derive 
analytical estimations for the matrix elements of the evolution operator in the long- 
time limit. In this way, approximations for the ground-state energy, as well as for 
the expectation of a generic operator in the ground state of a lattice system without 
sign problem, have been obtained as the solution of a simple scalar equation. The 
results presented in [H] were based on the application of a central limit theorem to the 
rescaled multiplicities of the values assumed by the potential and hopping energies in 
the configurations dynamically visited by the system (see also Ref. |S] for details on 
the corresponding probability density). 

The approximated scalar equation for the ground-state energy obtained by using 
the central limit theorem contains only the first two statistical moments of the 
multiplicities above mentioned and, as anticipated in [S], suggests that it is a second 
order truncation of an exact cumulant expansion. In this paper, after reviewing in 
detail the approach developed in we perform a large deviation analysis of the 
relevant probability density and obtain this exact cumulant expansion. In principle, 
if all the cumulants are known, we provide a scalar equation whose straightforward 
solution is the exact ground-state energy of the system. In practice, we measure the 
cumulants by a statistical sampling and only a finite number of them is available. 
The corresponding truncated equation provides ground-state energies whose level of 
approximation depends on the values of the Hamiltonian parameters and on the size 
of the system. However, the convergence as a function of the number of cumulants is 
rather fast and, as shown in some example cases, the use of the first few cumulants 
provides results indistinguishable from those obtained by exact numerical simulations. 

The spirit of our approach is different from that of Monte Carlo simulations. In 
the latter case, the accumulated statistical data provide information on the specific 
system under investigation and new simulations must be performed each time the 
value of a single parameter of the model is changed. On the other hand, our approach 
is semi-analytical in the sense that we have an exact expression for the ground-state 
energy, in which some statistical data, namely the cumulants of the multiplicities 
above mentioned, must be provided as an input. However, these multiplicities, and so 
the corresponding cumulants, do not depend on the values of the parameters of the 
model but only on the structure of the Hamiltonian operator (form of the hopping 
and interaction terms, size of the system, and so on). For instance, in the case 
of N spin-up and N spin-down particles with on-site interaction of strength 7 the 
multiplicities depend on the spectrum of the possible on-site pairs, 0, 1, 2, . . . , A^, not 
on 7 that represents, therefore, a parameter with respect to which our results are 
analytical. The analyticity with respect to the parameters of the Hamiltonian allows 
the determination of the expectation of a generic operator in the ground state of the 
system via the Hellman-Feynman theorem and suggests other potential applications 
that we will discuss at the end. 
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The paper is organized as follows. In Section El we briefly review the probabilistic 
representation of quantum dynamics for generic lattice systems. In this way, the 
ground-state energy is obtained as the expectation of a suitable stochastic functional. 
In Section 131 we summarize our main result, namely an exact formula for the ground- 
state energy, and provide an outline of its proof. In the four Subsections of the core 
Section 01 wc develop in full detail this proof via the analytical determination, in 
the limit of an infinitely long time, of the expectation of the stochastic functional 
providing the ground-state energy. In Section Wa\ we decompose the expectation in 
a series of canonical averages of weights, which are calculated in Section 14.21 The 
canonical averages are evaluated via a cumulant expansion theorem in Section [4. 31 and 
are finally resummed in Section 14.41 where the exact scalar equation for the ground- 
state energy is presented. In Section we discuss the numerical evaluation of the 
cumulants which appear in our formula for the ground-state energy. In the same 
Section, we also study some example cases and compare the results from our formula 
with those from exact numerical calculations. Finally, general features of our approach 
and future applications are summarized and discussed in Section 

Up to Section O we carry on the development both for hard-core bosons and 
fermions, whereas successively wc limit the method to hard-core bosons postponing 
the fcrmion case to a later work. 

2. Exact probabilistic representation of lattice dynamics 

In this Section we review the exact probabilistic representation of the imaginary- or 
real-time dynamics of a system of bosons or fermions on a lattice, see Ref. |3| for a 
detailed description. This representation is at the basis of our approach to study the 
ground-state properties of the system in a semi-analytical way We concentrate 
on a simple exclusion dynamics, in which multiple occupancies of the lattice sites are 
forbidden, i.e. fermions or hard-core bosons are considered. 

Let A be any finite lattice with cardinality |A| and 1,2,..., |A| some total ordering 
of the lattice sites. The dynamics of a system of hard-core bosons or fermions on this 
lattice is effectively described in terms of commuting or anticommuting destruction 
operators with the property (cj,j)^ = 0. Here, i = 1,2, ... , |A| is the site index 
and a =][ the spin, or an arbitrary internal index. The system can be represented in 
terms of Fock states, n e F, n = (77,11,7111, . . . , n|A||, »7|a|x)j where Ui^ is the lattice 
occupation number of the site i with spin a taking the values or 1. The number of 
Fock states is 



operator V is arbitrary, e.g. for the Hubbard model V = X^jgA 'Ti^It^iT^Ii^ii- ^'^^ 
kinetic operator we assume the quadratic form 




(1) 



where Nn is the number of particles with spin a. 

We assume the system to be described by the Hamiltonian 





(3) 
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with rjij„,9ija G R and rjij^ > 0. The case Oij^ 7^ is obtained in the presence of 
a magnetic field. In principle, our approach can be extended to more general kinetic 
operators. The essential feature to be noted is that in the Fock representation V is 
diagonal whereas K is off-diagonal. 

In order to study the ground-state properties of the Hamiltonian H, it is sufficient 
to evaluate the long-time behavior of X]n("'l^^^*l'"'o)- In fact, the ground-state energy 
is given by 

Eo = lim -atlogV(n|e-^'|no). (4) 

t^oo ^ — ^ 
n 

The quantum expectation of a generic operator O in the ground state of H can be 
obtained by evaluating the ground-state energy -Eo(0 of the modified Hamiltonian 
H -f and using the HcUman-Fcynman identity 

Let T(j be the set of system links with spin cr, i.e. the pairs («, j) with i ^ j and 
i, J G A such that ijija 7^ 0. To each link {i,j) with spin a we associate the value 

A,j,(n) = (n © 1,, © l,<,|e"'-4c^, + e*''-"cj,c,jn), (6) 

where \i„ = (0, . . . , 0, 1;^, 0, . . . , 0) and n ® n' = (n -I- n') mod 2. We may have 
lAijffl = 0, 1 only. A link (i, j) with spin a is called active if Ayo- 7^ 0. For 9ij„ = 0, in 
the case of hard-core bosons an active link can take only the positive value Xija ~ 1, 
whereas in the case of fermions an active link is cither 1 or —1, depending whether 
the number of particles present between the sites i and j, X]fc=i+i '^fco-, is even or odd. 

Let {iV*^-^} be a family of |r| = U Fj | independent Poisson process with jump 
rate p. These processes are in a one-to-one correspondence with the links of the lattice. 
At each jump of a Poisson process relating sites i and j with spin a and taking place 
at a given configuration n, a particle of spin a moves from site i to site j or vice versa 
if |A,jcr(n)| = 1, whereas the lattice configuration n remains unchanged if Xij^{n) = 0. 
By arranging the jumps according to the times s^, k = 1, . . . , iVj, at which they take 
place in the interval [0,i), we define a trajectory as the sequence of configurations 
ni,n2, . . . ,nNt generated from a chosen initial configuration tiq by exchanging, at 
each jump of the Poisson process N^^^, the occupations of the sites Uia- and nja-. The 
number of jumps Nt is, of course, a random integer associated to each trajectory. Let 
{ikTik) with spin Uk be the link jumping at the time Sk. By putting for brevity 

Afc = Ai,j,^,(nfe_i), (7) 

Vk^m^o^a^, (8) 

we define with Ai , A2 , . . . , A^Vt , and ryi , 772 , . . . , t^atj , the sequences of the corresponding 
link values and hopping parameters. To each trajectory, we also associate the 
sequences, Aq, Ai, . . . , A^^ and Vb, 14 ■ • ■ , Vat^, representing the number of active links 
and the potential energy of the visited configurations 

= E E \KMk)\ (9) 

Vk = {nk\V\nk). (10) 
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For later use, we also define 

Tk - AkVk+i/^, (11) 

where e is an arbitrary reference energy. 

At any finite time t, the following exact probabilistic representation holds jS] 

(n|e-^*|no) = E(<5„,„„^A1^J, (12) 

where the stochastic functional Mt is defined as 




.i{sk—Sk-i) I — ViVt(*-SiVt) 



(13) 



if iVt > and A^^„ = el'"l''*e~^"* if A^t = and the symbol E(-) means expectation 
over the |r| independent Poisson processes. In Eq. (|13|l . we put sq = 0. According to 
this representation, we have 

|no) = E(X:,J, (14) 



n e 



so that we can evaluate the ground-state energy as 

i?o= lim -ailogE(A^^J. (15) 

In the following, we will be able to reduce the study of the expectation of the 
stochastic functional (|13|l to the study of simpler canonical averages over the stochastic 
trajectories no,ni,n2, ■ ■ ■ in which ^ Wfe-i for any k > 1, i.e. trajectories in 
which only jumps corresponding to active links are considered, and the jumping times 
are disregarded (integrated out). These sequences of configurations no, rii, n2, . . . 
constitute a homogeneous Markov chain in a finite state space, F, with transition 
matrix P defined by 

{^(n)-i if 3 cr e {T,i} and {i,j) e T„: n' = n© 1^,, © 1^-^ 
and X^jaln) 7^ (16) 
otherwise, 

where A{n) = X]cr=ti, S(i j)er I'^uo-('T')| is the number of active links of the 
configuration n. Note that J2n' Pn.n' = 1 according to the fact that Pn^n' is the 
probability for the transition n n' . Finally, the canonical averages mentioned 
above are introduced for trajectories with a finite number of jumps in the following 
way. Given an initial state tiq and an application / : F^"*"^ — > C, function of ng and 
of the consecutive A^ configurations rii, n2 . . . , njv, we will indicate with {f) called 
canonical average, the average of / sampled with respect to the measure induced by 
the transition probability l|16() iterated A^ times. Similarly, we can consider canonical 
averages in which the initial state is not a single Fock state, no, but an ensemble with 
distribution ttq. An ensemble of particular interest is the invariant measure tt, defined 
as the left eigenvector of the transition matrix, -n^P = tt"^. It is simple to verify that 
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3. Main result and outline of the proof 

Evaluating the expectation E(A^^^) over the detailed realizations of the stochastic 
processes specified above can be done numerically by a Monte Carlo method j^j ■ The 
efficiency of the corresponding numerical algorithm is discussed in detail in Ref . \7\ . 
In this paper, we are interested to obtain analytical expressions of E (AlJ^^J, which 
are asymptotically exact in the limit t — > oo, and provide an exact formula for the 
ground-state energy by using Eq. l[T5|l . 

Let us call Jif, Y and ^ the sets of all the possible different values A, V and T 
assumed by Eqs. (0), and Hll() . respectively, along a trajectory no,ni,n2, ■ ■ ■ 
formed by infinitely many jumps. Let m^, my and he the corresponding 
cardinalities. Since any configuration can be obtained from any other one by a finite 
number of jumps, i.e. the Markov chain of the trajectories is irreducible, the elements 
in the sets Jif, 'V and 3" do not depend on the initial configuration np. Moreover, the 
value of their elements and, in particular, their number depend only on the structure 
of the Hamiltonian operator, not on the values of the Hamiltonian parameters. 

As we shall show, a crucial point is that, if we consider the conditional expectation 
in which the number of jumps N is fixed, E (A^^gl^Vt = A^), and integrate over all the 
possible jump times, what matters in determining the value of E (A^^„|A4 = 7V), is 
not the detailed sequences of the visited configurations but only the multiplicities, or 
numbers of occurrences, iV^, Ny and Nt of the variables A, V and T, respectively. 
Explicitly, for a trajectory with N jumps these multiplicities are defined as 

Af-l 

= E -^A.^A, (18a) 

iV 

^T=E%.T, (18c) 

with Afc, Vfc and Tk given by Eqs. 0, ^ and The expectation E (TW^jTVt = iV) 
is thus reduced to an average over the variables N\, Ny and Nt- In other words, 
after the stochastic times have been integrated out, the knowledge of all the canonical 
moments 

(7V„,...A„J^, (19) 

where e = 1^ U ,T U i = 1, . . . , fc, determines completely the expectation 
£{Ml^Nt=N). 

Let us consider the case of hard-core bosons in the absence of magnetic fields, i.e. 
the case for which ^ = {1}. As we will prove, after the integration of the stochastic 
times, the conditional expectation E (A4^,jA^t = A^) becomes 

E {Ml^N, ^N) = lwN{t) n ) , (20) 

where WAr(i), called weight, is a function that depends only on the multiplicities of 
the potential (|18fe|l , whereas the other factor is a purely kinetic function that depends 
only on multiplicities H18c|l . 
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In order to understand the behavior of Eq. H2()|l . and then of Eq. (|14|l . let us 
start to consider a non interacting system with V = Q. In this case we have 1^ = {0} 

and the weights have the simple exact expression W^^(t) = e^t^ /N\. Therefore for 
E {M*n^\Nt = iV) we are left with a residual canonical average over the variables Nt 
and we get 

E {Ml^N, -N)=^^ (ei:^'-(-)--)^ ^ (21) 

where a is a suitable value coming from the integration over Nt- 

The result in Eq. ^ can be immediately obtained by the rough estimate 
(exp[^j, log(T)A'^r])jy ^ const^, based on the boimds exp[A^log(Tniin)] < 
(exp[^j, log(T)A'^r])jv — cxp[iVlog(ri„ax)], which in turn follow from the 
normalization X]Te5' ~ N. By using a cumulant expansion theorem [S], however, 
it can be shown that this result becomes exact for ^ cx) if we assume the existence 
of the following rescaled cumulants of the variables Nt 

^^Z..,n^^^l^j^{Nn---Nn)^N\ (22) 

where with we indicate the cumulants, or connected correlation functions, 

sampled with respect to the measure induced by the transition matrix (|16|l . The 
connection between these cumulants and the statistical moments introduced above, 

is well known. For fc = 1 we have (A^tJ^^^ = {Nt,)n, for k ^ 2 (^TiA^Ta)!^^ = 
{Nt,Nt^)j^- {Nt,)^ {Nt^)j^, and so on. 

Since the Markov chain with transition matrix H16|l is finite and irreducible, for 
k = 1 the existence of the limit (|22|l is assured by an ergodic theorem [0]. On the other 
hand, for fc > 1, the existence and the finiteness of the rescaled cumulants (|22|l . or of 
the more general limits (|26|l , follows from the exponential decay of the correlations of 
local functions of the configurations along the Markov chain. For the T variables, this 
property amounts to say that for any k do exist positive constants bk and Nc{k) such 
that for any N 



{^Xt, (rih, )---Xt, < bk exp j 



(23) 



where Xr('^) ~ ^T(n),T is the characteristic function of the states T and /imax and 
^min are the maximum and the minimum of the indices 0<hi<N,l = l,...,k along 
the chain. A numerical check of Eq. (|23|l for k ~ 2 will be shown in Sectional 

Since the rightmost expression in Ea. (|21|l has, as a function of N , an exponentially 
leading maximum at iV = eat, we see that for t large an important consequence follows. 
The full expectation E (Al^^) = Sw'=o ^ {-^nj^t ~ takes exponentially leading 
contributions from terms with ~ eat. Therefore, for t cx) we have that i) 
the saddle point technique used to evaluate the canonical averages over iV^ becomes 
exact and ii) sl further saddle point integration can be used to exactly resum the full 
expectation E (A^^^J. In conclusion, we obtain an exact expression for the ground- 
state energy, namely hmt^oo ~dt log E (M^J. 

For V 0, the above described scenario remains essentially unchanged. In fact, 
even though in this case the integration of the stochastic times cannot be done exactly, 
it can be performed by a saddle point approximation, which becomes asymptotically 
exact for TV — > oo. Independently of their exact value, as we will prove later, the 
weights are bounded by wj^' (<)e-^"""'* < WN{t) < >vj^^(t)e-^"""* and behave, as a 
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function of TV, similarly to the non interacting case. The conditional expectations H2()|l . 
which now reduce to residual canonical averages over both the variables Nt and Ny, 
can be evaluated analogously to Eq. (|21(l . The result again implies that the saddle 
point integrations we perform to evaluate i) the weights, ii) the residual canonical 
averages, and Hi) the sum of the series X]^=o ^ (-^tioI-^* ~ become exact in 

the limit t ^ oo. The conclusion is that for hard-core boson systems in the absence 
of magnetic fields, i.e. = {!}, the ground-state energy E^b is determined by the 
following scalar equation involving all the cumulants of the Markov dynamics 

oo 

Efc! E E ..,a."o,(£;OB)...«a.(i?OB)=0, (24) 

fc=l ' axe.^C Q&e.-^ 

where = "^0=5^, and the vector = (... uy •■ •) with V <^ 'Y and 
T e 5^ is defined as 



u^{E^b) = (•••- \og{{-E^B + y)le\ logT .. .). (25) 
The asymptotic rescaled cumulants Eai,...,afc are defined as 

(26) 

where with (•)^^ we mean the cumulants, or connected correlation functions, of the 
multiplicities A^q, a S , sampled with respect to the measure induced at N jumps 
by the Markov chain with transition matrix 1)16(1 . 

For F = 0, Eq. H24|) can be solved explicitly and we get the following exact 
expression for the ground-state energy iJg^^ of a non interacting hard-core boson 
system 



i:.(0) 

^QB = -eexp 



E ••• E sS....T.iogri...iogT. 



^ k. 

,k=i Tie,3r TfcG^^ 



(27) 



For a generic V, Eq. (|24|l . with the series truncated at an arbitrary finite order 
^max, can be solved numerically in a straightforward way by using the bounds 
-E'oB — ^OB < 0. 

Equation (|24(l allows, via the Hellman-Feynman theorem also to evaluate the 
expectation of any other operator in the ground state of the chosen Hamiltonian (see 
Ref. |S] for more details). Moreover, it is clear that the cumulants S^*"'^ depend only 
on the structure of the Hamiltonian operator, not on the values of the Hamiltonian 
parameters. Therefore, once the S'^'^^ are known, all the evaluated ground-state 
expectations are analytical functions of the Hamiltonian parameters. 

Our analytical formula for the ground-state energy ()24|l rests on the knowledge 
of the asymptotic rescaled cumulants. In some special cases, it could be possible 
to evaluate, at least approximatively, also these cumulants analytically. We refer to 
Section [S] for a detailed discussion on how determining the cumulants in a numerical 
way. In that Section, some example cases arc worked out explicitly. Finally, we refer 
to Section El for a scenario of possible applications that exploit the analytic character 
of our approach. 
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^■l. Canonical decomposition of the expectation 

To evaluate the expectation E(A^^^J, we decompose it as a series of conditional 
expectations with a fixed number of jumps (canonical averages) 

oo 

^{Ml^) ^Y.^^^'r.^Nt=N). (28) 

In the canonical averages the stochastic jump times arc easily integrated out. In 
fact, multiplying the stochastic functional (|13|l with the condition Nt = N by the 
infinitesimal probability, 

dPN = e-^'^^'^pdsi e-|ri(^=-^iVrfs2 . . . e-l^l(''"-""-iVf^sw e-I^K*-''"^ 

to have the jumps 1,2, . . . ,iV of the independent Poisson processes in the intervals 
[si, si + dsi), [s2, S2 + ds2), • ■ • , [sni sn + dsj\j), respectively, and integrating over the 
times si, S2, . . . , s^v, we get 

E {MljN, = iV) = E 5(;)/cWw«(0, (29) 

where fijv = QnItiq) is the set of all possible trajectories with N jumps branching 
from the initial configuration no and 

S^'^ = ^ ^2 ' • • ■ (30) 
/C(;)==e-V%«...,y«, (31) 

Equations (|3l)|l . H31|) and H32|l define three dimensionless quantities related, 
respectively, to the sequences of the signs (more generally, of the phases), of the 
hopping parameters and of the potential values associated to the rth trajectory with 
iV > jumps. We have 5^'"' = = 1 and w''^\t) = e'^"* in the case N = 0. The 
quantities W^\t) are positive definite and will be called weights in the following. 

From Eq. (|3()|l . it is clear that only the trajectories formed by a sequence of active 
links contribute to the sum in Eq. 1)29(1 . Hereafter, therefore, we restrict ^In to be the 
set of these eff'ective trajectories with N jumps branching from no and we exclude the 
value A = from the set 

The sum over the set can be expressed as an average, (•)^, over the trajectories 
with N jumps generated by extracting with uniform probability one of the active links 
available at the configurations ng, ni, . . . , n7v-i. The probability associated to the rth 
trajectory generated in this way is pj^^ = nl-io^ ^/^^k^ ^^'^ have 



N-l 

v^Qn rG^N k—Q 



= /sNWNit) YlnS , (33) 

\ fc=0 / N 
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where (•)^ = X^rsnjv ' P'n ■ Note that X^rsOjv P^n ^ 1- -^y using the definitions (|18Qjl 
and H18c|l . we rewrite the canonical averages as 



E {Mi^}N, ^N)^ lwN{t) n n ^""n 



(34) 



N 



In the foUowing, we will consider only hard-core bosons in the absence of magnetic 
fields. In this case, Jtf = {1} and we are left with averages 
of positive definite quantities 

E {Ml,\N, ^N) = /w^{t) l[T^A . (35) 
4-2. Evaluation of the weights 

In this Section, first we find a recursive differential equation for the weights >VAr(t). 
Then, taking the Laplace transform of this equation, we realize that the weights 
depend only on the multiplicities Ny of the potential, not on the detailed sequence 
Vb, Vi, . . . , Vat- Finally, we use a saddle-point technique on the complex plane to obtain 
an explicit expression for the weights, which becomes exact in the limit N — > oo. 

Equation H32() . which defines the weights >VAr(t) for > 1, can be rewritten as 
(for simplicity, we omit the trajectory index (r) but we stress the dependence on the 
sequence of the potential values) 

WNit; Vo, Vi,..., Vn) = GNit; Vo,Vi,..., Viv)e-^~*, 

where 



GN{t;Vo,Vi,...,VN)^e^ dsi / ds2 . . . / dsN9{s2 - si)9{s3 - S2) . . . 

Jo Jo Jo 

X 9{sN - sjv-i) exp (Aisi + A2S2 + . . . + Ansn) , 

with Afc = Vfc — Vk-i, k = 1,2, . . . , N . By evaluating the derivative of G'Ar(t), for 
> 0, with respect to t, 

dtGN{t; Vq,Vi,..., Vn) = eGN-i{t; Vq,Vi,..., l/^-i) cxp (A^i) , 

where Go(t; Vb) = 1, we find the following recursive ordinary differential equation for 
Wiv(t) 

dtWN{t; Vo,Vi,..., Vn) = eWN-i{t; Vo,Vi,..., Vn-i) 

~VNWN[t;Vo,Vi,...,VN). (36) 

Since Eq. 136|l is linear in yVAr(i) and WN-iit), it is convenient to introduce the 
Laplace transform 



yVN{z) = / dte-'^'WNit), 2 e C. 

By observing that WNi^) = for > 0, Eq. (|36|1 reduces to the following recursive 
algebraic equation for Wat (2) 

zWN{z)^eWN-i{z)-VNWN{.z), (37) 

from which we get 

WN{z) = e{z + VN)-^WN-i{z). (38) 
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We recall that for TV = we have >Vo(i; Vq) = e"^"* and, therefore, Wt){z) = [z+Vq)-^. 
By iterating Eq. ()38|l . we arrive at the solution 

k—v 



which, in terms of the multiphcities (|18&|1 . takes the form 

WN{z)^e^ \{{z + V)-''-. (39) 

This expression shows that, for any value of N , the weights depend only on the 
multiplicities Ny, i.e. WN{t; Fq, Fi, . . . , Vn) = WN{t; {Ny})- 

The explicit inversion of the Laplace transform (|39|l can be done analytically only 
for N large. However, as we will immediately prove, this is the limit we are interested 
in. In fact, in the case ee 0, we have WN{t) W^^(0 = e^t^ /N\, which is easily 
obtained by direct integration of Eq. %VZ\ . For t large yV^''(t) has, as a function of 
A'^, a maximum at N{t) = et and, around this maximum, is well approximated by 
a Gaussian of width y/N{t). For V generic, let us indicate with Vnun and Knax the 
minimum and maximum elements of Y such that Ny 7^ for the trajectory we are 
considering. In Eq. (p?^ we thus have T4iin < Vk < Vmax, fc = 0, 1, . . . , and 

We conclude that the weights arc bounded by 

Wj,°)(t)e-^— * < WN{t) < Wj,°)(i)e-^™*, (41) 

so that, for t — > 00, exponentially leading contributions are obtained from values of N 
in the range [et — ^/ei, et + \/ei] . 

According to Eq. WnIz) has ttin < my poles at the points Zy = —V such 
that V G y and Ny ^ 0. For A^ sufficiently large, the number tojv of these poles 
approaches m^/. For any finite value of A^, the poles of WAr(z) arc confined in the real 
segment [— Knax, — ^lin]- Recalling the rule for the Laplace inverse transformation, 
we have 

Wjv(t) = 7^ f dze^'WNiz), 

where the integration contour F can be any line parallel to the imaginary axis and 
contained in the analyticity domain of the Laplace transform. In our case, F must 
be in the domain Rez > — Vmin- By Jordan's lemma, the contour can be closed to 
infinity in the left half-plane Rez < —Vmin without changing the integration result. 
Finally, by Cauchy's theorem, F can be deformed into any other anti-clockwise closed 
contour F' still containing all the poles^zy. 

By using the expression (|39|l for yVAr(z), we write its antitransform as 

Wjv(t) f dz cxp [A^^(z)] , (42) 

Zme Jy 

where 

zt s:^ Ny Z + V 

For t large and N ^ t, since also Ny ~ A^ due to the ergodicity of the trajectories, 
we can evaluate the complex integral 1421) by a saddle-point technique. Let us call zq 
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the point of the complex plane where ^p{z) is stationary, i.e. Dz^p^zq) = 0. We deform 
the contour T into a new one F' (steepest descent) such that, for z e F', Reip^z) has 
a maximum at zp, whereas Imi^(z) is constant, at least at the first order in its Taylor 
expansion around zq. Provided that the contour F' remains in the analyticity domain 
of ip{z), the following standard result applies [TU| 



^N{t)^-^J—j^^^^exp[ia + N^{zo)], (44) 
27riey N\diip[zo)\ 

where a is defined as 

a=|-^Arga?^(zo) (45) 

and represents the angle formed, at the saddle point zq, by the new contour F' with 
respect to the original one F. 

Now we calculate the first and second derivatives of </?(z) and the saddle point zq. 
From Eq. gSJ we have 



In terms of a;o = Rezg and j/o = Ihizq, the saddle point equation, dz^p{zf)) = 0, reads 

t ■r-^ Ny Xq 

As Ny > 0, Eq. Ij486|l is satisfied only by j/o = 0. Hence, we are left with the following 
equation for xq 

E = t- (49) 

For any i > 0, Eq. H49|) has tum solutions. The first m^r — 1 solutions, ordered 
according to their increasing value, are in the range — VJnax < a;o < — Knin, whereas 
the last one is such that xq > — Knin- This is the only solution compatible with the 
condition that the integration contour F, passing through zq = xq + iyoi is contained 
in analyticity region of Lp{z). Therefore, for any t > 0, we have one and only one 
saddle point determined by Eq. with the constraint .tq > — Knin- 
The fact that yp = also implies 



so that we have 

ArgaXzo)=0, (51) 

i.e. a ~ 7r/2. The integration contour F, therefore, has to be deformed into a 
new one F' parallel to the real axis. At a first sight, this kind of deformation is 
incompatible with the condition that F' strictly contains all the poles Zy located on 
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Figure 1. Integration contours F and F' in the complex plane z = x + iy, which 
are used to evaluate the Laplace antitransform of WAr(2). The saddle point xq 
and the poles of Wn{z) are indicated by open and filled circles, respectively. 



the real axis. However, the largest solution of Eq. H49|) is bounded from below by 
Xq > —Klin + Ny^i^/t. This means that for i ^ oo with Ny^^^ ^ N ^ t, the distance 
between the saddle point and the closest pole at z = — V^nin is finite. Therefore, we 
can take F' parallel to the real axis only in a finite neighborhood of the saddle point 
as shown in Fig. ^ In fact, the contribution to the integral due to this neighborhood 
of Zq is exponentially leading in N with respect to the rest of the contour. 

In conclusion, for t oo with N ^ t, the integration on the contour F' chosen as 
described above becomes asymptotically exact and, in this limit, we have the following 
exact expression for the weights 

gOOot-Y: ve r \o^[(xq+V) / e] 
V ^t, Xo> -Klin. (526) 

^-^ xo + V 

This expression is semi-analytic in the sense that the simple Eq. (|52b(l , which provides 
Xo to be inserted into Eq. 1)52 all , must be solved, in general, numerically or recursively. 

For 1^ = 0, we can verify that the above expression for the weights reproduces the 
exact result W^^\t) = e^t^/N\. In this case r = {0} and the solution of Eq. (|52b|l 
is 

N + l 

Xq = . (53) 

Inserting this value into Eq. (|52Qjl . we have 

^(0) ^ exp[jV + l-(jV + l)log[(jV + l)/ef]] 
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1 exp(7V + l) 



e"t". (54) 



V2^(Af + 1)^+1/2 

By recalling Stirling's formula 

N\ ~ V2^iV^+i/2e-^, (55) 

which derives from a saddle-point evaluation of the Gamma function as well, we see 
that Eq. (|54|l is just the Stirling approximation to e^t^ /N\. 

4-3. Canonical averages via large deviation analysis 

To evaluate the canonical averages it is useful to introduce the frequencies, Vy = 
Ny/N, V € and i/j, = Nt/N, T € ^, which for N large become continuously 
distributed in the range [0, 1] with the constraints 

Note that, for N large, we do not distinguish the different normalizations. A'' + 1 and 
N ^ of Nv and Nt: respectively. When possible, we will use a compact notation in 
terms of the vectors jj, and i/, which have m = m-y + m,^- components indicated by 
a Greek index a S J/iC = Y U .!7 and are defined as = (■ • ■ Nt ■ ■ ■) and 

= {. . . Vy . . . \ . . . Vj, . . .) . We have 

/X = Nu. (57) 

For later use, we also define ~ (• • • ^ log[(2;o + ^)/e] logT .. .), = 

(. . . (xo + V)~^ .. .) and = (. . . (zq + 1^)"^ .. .)• Note that the 

vectors M, V and w depend on u through = xq{u) and v = —OxgU, w — —dx^v. 
Finally, we will take advantage of a scalar product notation. For instance, we rewrite 
the saddle-point Eq. (|52b(l as {v,v) = t/N. 

By using the result given by Eqs. H52ajl and 1)52 we express the r.h.s. of Eq. (|35)1 
in the following explicit form 

\ pXot+{fi,u) 

The probability VNifJ-) is given by the fraction of trajectories branching from the initial 
configuration tiq and having after N jumps multiplicities /x. According to Poisson's 
summation formula 

E/(/^)= E / rfMe'('=^'^V(M), (59) 

the sum over fi in Eq. (|58|l can be transformed into a series for k e Z™ of integrals 
over the same variable in the presence of the oscillating factors exp[i(fc, /x)]. As we 
will check at the end of the calculation, in the limit of t large all the terms fc 7^ of 
this series are exponentially damped with respect to the term fc = 0. In this limit, 
therefore, we will not distinguish the sum over in Eq. H58() with the corresponding 
integral. 

In we have evaluated Eq. 1)58)1 for t large by using a central limit theorem 
for Markov chains. Although this theorem applies rigorously to each rescaled sum 
Na/^/N = Va\/N, a € M', it provides an approximation when, as in Eq. 158)) . 
variables like v^N are involved. As anticipated in |^, the integrand in Eq. 1)58)) is 
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sensitive to the large deviations of 7'Ar(/i.) from the central limit behavior and to 
obtain more accurate estimates we need to consider a development in terms of the 
associated cumulants (connected correlation functions) |H]- 

Before proceeding with this analysis, we observe that the constraints H56|l give 

important summation rules for the cumulants. Let us indicate with {vai ■ ■ ■I'Qfcliv^ 
ai, . . . , Qffc G J^, a cumulant of order k. In Appendix |B| we demonstrate that 

for fc = 1, whereas for fc > 1 

E (''"I •■•^"Jw = E ('^ai-.-i^ajJ^^ =0. (61) 

These rules provide a basic test for the statistical measurement of the cumulants 
themselves. A sampling that aims at measuring the cumulants with a given accuracy, 
will have to satisfy Eqs. (|60|l and H61|) with the same accuracy. 

As customary in the framework of a large deviation analysis, we are interested to 
get information about the density Vn{Ni>) in the limit of N large and v finite. By 
introducing the Fourier anti-transformation 

VNipi) = j dqe'°s[^"('j)]-^('^.'^), (62) 

Vn{(i) being the Fourier transform of 'Pn{iJ'), Eq- (|58|) becomes 

(w^it) ^T^-^ = (^^) " Jd'^j dqe^'^^-^^'^^Riu), (63) 

where ^(f, q) takes into account the exponential behavior of the integrand, whereas 
i?(i/') is a smooth function. For brevity, we omit the parametric dependence of (f) and 
R on t and N. Explicitly, we have 

c^iu, q)=x,^ + {v,u)-i{v,q) + l^lM^ (64) 

R{u) = ^ (65) 

As well known, the cumulants arc related to logPjv(q) through the relation 
log^Ar(q)= log(e'('^''')' 



N 

1 (^) 

= Em((/^'^'?)'>J' 

k=l 

(<--)'>:'■ w 

k=l 

Note that for any given N ^ due to the inequalities {^ai ■ ■ ■ IJ-Ok) n !i A^'^, valid for any 
fc, and due to the asymmetry 7'Ar(/x) ^ 7'Ar(— /x), the series in Eq. H66(l converge for 
every q £ C™ (see, for example |H]). 

We evaluate the integrals in Eq. H63|l by a saddle-point calculation (actually 
Laplace's method) with respect to the variable (i>',q) £ R'" x R"\ The derivatives of 
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4>{i', q) with respect to qa and u^, a ^ are 

where we have used {i/, v) = t/N. Therefore, the stationarity condition for q) 
imphes 

^ N'-' I, . ,,-1 \W 



which, in turn, reduces to the following equation for the saddle-point frequencies i/^p 

°° /V*-''! / \ (c) 

'^"'' = E(fc3I)!((-'"(-o('^")))'"-">^ ■ (67) 

Hereafter, we will add the superscript '^p to a function of u to indicate the value 
of this function for v = f'^P. By approximating the function (/)(i>',q) with its second 
order Taylor expansion around the saddle point {i^^^, q^^) and performing the resulting 
Gaussian integral with respect to the variable {u — v^^, q — g'^P), Eq. becomes 



Wiv(0 YlT^'A = Cn exp [V'(iV)], (68) 
Te.sr I ^ 



where 



V'(iV) = A^0(z^"P,q"P 

A:! X^" ' Wat 



fc=i 

= , , (70) 

v/|detV2^(iv^P,q^P)| 

V^(?i>(i>', q) being the Jacobian of (f) with elements 

fe=2 ^ ^' 

d^^d^fj(t>[u,q) = 



dq„d^f34>{v,q) = d^^dqij(j){v,q) = -iSa,0, 
for a, /? G It is convenient to introduce two mxm matrices S and A with elements 

So,;3= -a,„a,,</>(i.^p,0 (71) 

^o,;3= -a,„a,„0(/y^P,O. (72) 
In terms of these matrices Eq. (|70ll becomes 

Cn = ^ = , ^ . (73) 

v/|det(l + SA)| ^27riVe2(i^sp^ ^^^^p) 

Note that, in the case V = 0, the matrix A is uniform and, due to the summation rules 
H6U|I and H61|) . HA = 0. In general, the same summation rules imply det(SA) = 0, 
so that dct(l + SA) ~ 1 up to second order terms in A. 
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4- 4- Resumming the canonical series 

In this section, we evaluate the expectation E (AlJ^^) by resumming the scries of 
Eq. (|2H1)- We will substitute the series over N with the integral 



(74) 



and we will perform the integration again with a saddle-point technique. At the end, 
we will demonstrate the asymptotic exactness, in the limit t oo, of the saddle-point 
integration, as well as of the substitution of the series with the corresponding integral, 
by showing that in this limit Eq. (|74|l takes exponentially leading contributions from 
values of N in the range [et — Vet, et + \/ei] . 

Before evaluating the stationarity condition for 'ijj(N), an important comment is 
in order. For any finite value of N, let us introduce the rescaled cumulants of order k, 
in a compact notation , which is defined as the tensor of rank k with components 

Sl^^'^L. =7V^-V;.„, ...i^„J^\ (75) 



ai,...,afc € J^f. Let S^*^-* be the tensor of the asymptotic values of the rescaled 
cumulants in the limit N oo 

E« = hm 7V'^-i(^.„, = hm -'^r m \W 



. V, 



aki N 



N 



(76) 



These limits exist and are finite since the irreducible and finite Markov chain formed 
by the evolving configurations has a finite correlation length Nc with respect to the 
number of jumps (see Eq. ^^). Up to corrections exponentially small in N/Nc, for 
N ^ Nc we can use the effective rule 



0. 



(77) 



By using the rule H77(l . the derivative of iIj{N) with respect to N is 



fc=i 



N 



-^^(fe-1) 



N 



On the other hand, from Eq. (|67|l we have 



so that, by using N{iy'^'^,v^'^) ~ t, we find 



fc=l 



fc! 



(c) 
N 

(c) 
N 



In conclusion, the stationarity condition, dNip{N^^) = 0, gives 



E 

k=l 



k-l 



kl 



N 



-0, 



which in terms of the rescaled cumulants (|75|l reads 



fc=l 



= 0. 



(78) 



(79) 



(80) 
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Equation (|80|l determines the saddle point N^^ for the chosen value of the time t. 
From Eq. (f7H|l and using the rule (|77jl . for the second derivative of ip{N) we find 

92^VW = -(i-^P,^^P)5rf, (81) 

where 

dj^xl^ = i-- ^ ' J -. (82) 

Equation H82|l has been obtained by evaluating the total derivate of the saddle-point 
equation N(y^P ^v^p) = t with respect to N , 



(iv^P, v^P) + N (a.^=piv^P, t;'P) - (j^'^P, ii^'^P) 
and by observing that, according to Eq. H67|l . we have 



0, 



°° A[k~l I , (c) 



(a..p.-,.-)= - ^ --- ^(.,«-)- i^^,v^r)^ 

k=2 ^ '' 

= - (Si;'^P,t;'^P). 

To evaluate the integral we approximate i}j{N) ~ i}j{N''p) + ia|^-0(A^^P)(7V- 
N^pY and Cjv — Cn^p- The remaining Gaussian integration gives the following result 



1 + tr(SA) 



|det(l + SA)| e(i^^P,t)'^P) 



(83) 



The matrices S and A arc defined by Eqs. (|71|l and H72() and their components 
explicitly read 

oo 

+ E ■•■ E (84) 

k=l aiG^ a^G^ 
;3 = , ° ^ , ■ (85) 

Equation (|83|) represents our final expression for the matrix elements (|14ll . It is 
valid for t large but finite with N^P{t) determined by Eq. H80() . It is simple to show that 
in the limit t ^ oo, Eq. H80(l may admit a solution only if N^P{t) grows proportionally 
to t. In fact, from Eq. I|526(l for arbitrary t and N wc have that 

sn ^ + 1 sn N +1 , ^ 

^K..„ < -^-o' + Knin < , (86) 

where Vrnin is the smallest value oi V G Y such that I'y > 0. Note that this value 
of V exists due to the properties < z/^p < 1 for any a G and J2a<£ r'^a' = 1- 
Therefore, for N = Nsp{t) all the components uy = — log[(a;g'^ + V)/e] diverge when 
t ^ oo if limt^oo N^^{t)/t = oo, whereas if limt^oo N^^{t)/t = the only diverging 
component is Uy . . In both cases, Eq. I|8()(l does not have solution for t —^ oo. In this 
limit, therefore, we must have N^P(t) ~ t. This condition, which from a physical point 
of view simply expresses the proportionality between the time and the length of the 
trajectories most contributing to the expectation, implies two important consequences. 
First, since the peak of the Gaussian at = N^P{t) moves to infinity linearly with t, 
whereas its width increases only as ^/i, see Eqs. (|81|l and H82(l . the result H83|l becomes 
asymptotically exact for t oo. Second, in Eq. H8U|I we can substitute S^^-*^) with 
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their asymptotic values S*^*^). In this way, Eq. (|8()|l can be read as a time-independent 
equation that determines the quantity 

Eqb = - lim XgPljY^jv.p.j) (87) 

as a function of the asymptotic rescaled cumulants S'*^'. A: > 1. By comparing 
Eqs. (|15|l and H83|) . we realize that this quantity is the ground-state energy of the 
considered system. In conclusion, for systems of hard-core bosons the ground-state 
energy Eqb is determined by the scalar equation 

oo 

Efc! E E 4';l.,a."o.(^^0B)...^io,(i?0i3)=0, 

fe=l ' aiG.>f ttfeS,^ 

where 

u^{Eqb) = (•••- log[(--BoB + V)le] logT .. .). (89) 

If we truncate the series in Eq. to the second order, we recover the Gaussian 
approximation of Ref. |5] (there S'^^-' and S*-^-* were indicated as V and S, respectively). 
Furthermore, as in the V = Q case can be solved explicitly and we get the following 
exact solution for the ground-state energy, E^g , of a non interacting hard-core boson 
system 



Kb = -eexp 



E^ E •■• E sS...,T.iogr,...iogT, 



k. 

= 1 Ti&Sr TfcG3^ 



(90) 



For a generic V , Eq. (|88(l . with the series truncated at an arbitrary finite order fcmax, 
can be solved numerically by the bisection method using the bounds E^^g < Eqb < 0. 

Equation (|88|) allows, via the Hellman-Feynman theorem (O, also to evaluate the 
expectation of any other operator in the ground state of the chosen Hamiltonian (see 
Ref. |Sj for more details). Moreover, it is clear that the cumulants S^*^^ depend only 
on the structure of the system Hamiltonian, not on the values of the Hamiltonian 
parameters. Therefore, once the S'*^^ are known, all the evaluated ground-state 
expectations are analytical functions of the Hamiltonian parameters. 



5. Numerical evaluation of the cumulants and example cases 

In this section, we discuss the numerical evaluation of the cumulants. We also apply 
our method to some example cases and compare the ground-state energies obtained 
by Eq. (|88|l with those from exact numerical calculations. 

In our approach, the starting point is the evaluation of the asymptotic values S'*^^ 
of the rescaled cumulants S'^^'*''\ According to Eq. ljA.7p . the latter are obtained 
in terms of the statistical moments {Na-^ . . . Na^) j^, which are easily sampled by 
generating AI random trajectories branching from the initial configuration no, i.e. 

M 

-mE^^^---^^.^ (91) 

(v) 

where Na is the multiplicity N^, a g Jif, of the pth trajectory. The number M of 
trajectories must be chosen larger than a critical value M^{N,k), which depends on 
the statistical precision e required in the evaluation of the rescaled cumulants Ti^'^'-^K 
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Figure 2. Absolute value of the normalized correlation function C'v =-y .v —-y {h) 
as a function of the jump interval h = h2 — hi for two different FNU hard-core 
boson Hubbard models. The straight dashed lines are the results of a fit in the 
range 10 < /i < 29 with a function of the form const X exp{—h/Nc). We have 
A^c = 4.80 and Nc = 6.14 in the 2 x 3 Ap = 3 and 4 X 4 Ap = 8 cases, respectively. 
Each correlation function has been evaluated by using A-I = 0.5 X 10® independent 
trajectories. 

The length N of the Af trajectories is chosen sufficiently large for the asymptotic 
behavior of S'^^''^^ to be established. This may represent a problem since the 
fluctuations in Eq. H91|l grow as A^'' so that M^{N, k) becomes huge for fc > 1. However, 
the evolving trajectories form a Markov chain having, see later, a finite correlation 
length Nc- Therefore, S^^'*^^ converges exponentially to S^'^^ with a characteristic 
length of the order of (fc + l)iVc, fciV^ if the initial configuration no is taken randomly 
distributed according to the invariant measure of the Markov chain. As shown in the 
example case of Fig. El for a large class of models we have observed that the correlation 
length Nc exists and grows no more than linearly with the size of the system. The 
evaluation of the cumulants is thus feasible even for large systems. 

A possible numerical limitation in the evaluation of the cumulants of large order 
is collecting and updating all the components of S^*^', whose number is m'' with to, 
the cardinality of Jff, that grows with the size of the system. Basically, this kind of 
difficulty is related to the CPU capability of simultaneously updating many addresses 
of memory, not to the CPU speed, and can be reduced by vectorization tools. The 
difficulty can be lowered also exploiting the invariancc of the components T.^a^^,,,,a^. 
under permutation of any pair of the a indices. In this way, only the components 
sS...,afc with cti < a2 < ■ . ■ < ctk arc sampled according to Eq. H91(l . This introduces 
an error in the summation rules (|60|l and H61|) . which, on the other hand, are identically 
satisfied if all the components of S^'^^ are sampled. However, more than a drawback 
this error represents an advantage, which allows to set the critical number Af^ of 
sampling trajectories in a simple way. In fact, the determination of the cumulants 
with a statistical precision e implies the summation rules (|60|l and (|61|l to be satisfied 
with the same precision. 

In the following, we report on simulations performed in the case of the first- 
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neighbor uniform (FNU) hard-core boson Hubbard model defined by the Hamiltonian 

(ij)ercr=n ieA 

where F = {(i, j) : i < j G A and i, j first neighbors}. We set the reference energy e 
to the value of the hopping parameter 77. We have considered two-dimensional systems 
with periodic boundary conditions, x Ly sites and N-\ = = Np particles per 
spin. For this model, the set Y has elements 0, 7, 2j, . . . , Npj, whereas the set is the 
collection of the possible values of the number of active links, e.g., ^ = {12, 16, 20} in 
the case of a 2 x 3 system with Np — 3 or ~ {8, 10, 12, 14, 16} for the same system 
with iVp = 2. 

To check the correlation properties mentioned before, we studied the connected 
correlation functions of order 2 

Ca.Ah2 - h,) - (x„ {n^,) Xp KJ)!;^ , (93) 
where > hi and Xa i^^) the characteristic function taking the value 1, if V{n) = a 
and a € 1^, or T{n) = a and a € =5^, and otherwise. The averages can 
be sampled as indicated in Eq. H91() by generating M independent trajectories with 
configurations n^i^\ h = 0, 1, . . . , iV and p ^ 1, . • . , M. For N sufficiently large, the 
correlation functions H93() do no longer depend on hi and /12 but only on the jump 
interval h = h2 — hi. In Fig. El we show the behavior of the correlation function 
(normalized to 1 at /i = 0) obtained by choosing a and /? equal to the potential 
value y = 7 for two different FNU hard-core boson Hubbard models. After an initial 
transient, Cv=j,v=-y{h) decreases as exp{—h/Nc). The measurement of the correlation 
length by a fitting procedure shows that Nc increases slowly with the size of the system. 
Similar results are obtained for different choices of a and /3. 

In Fig. 13 we show the behavior of Eob as a function of the interaction strength 
7 in a 2 X 3 lattice with Np = 2 and Np = 3, whereas in Fig. ^ we consider a 4 x 4 
lattice with Np = 5 and Np — 8. In these figures, we compare the energies obtained 
from Eq. by truncating the cumulant expansion at the order fcmax = 1,2,3,4 with 
the results from exact numerical diagonalizations (Fig. Ojl and quantum Monte Carlo 
simulations (Fig. EJ. For fcmax = 2, wc recover the results of Ref. 0. As expected, we 
obtain better and better agreement with the exact energies as the truncation order 
fcmax is increased. 

The number of cumulants needed to obtain a given approximation grows as the 
the interaction strength 7 or the lattice size |A| are increased. As explained in the 
previous Sections and anticipated in , this behavior is due to the form of the function 
f{fj,) to be averaged in Eq. (jSHJ- In fact, /(/x) involves multiplicities coupled with 
the potential values and with the number of active links, which, in turn, are related to 
7 and |A|. For increasing values of 7 and |A|, f{fi) becomes more and more sensible 
to the large deviations of the probability density PAr(/x) and cumulants of higher and 
higher order must be kept in the calculation. 

In Fig. which is an enlargement of Fig. 13 (left panel) and Fig. 01 (right panel) 
in the small 7 region of the systems at half filling considered there, we can better 
appreciate the convergence of the the solutions of Eq. for increasing values of 

fcmax, toward the exact energies. In all Figs. I3I5I the statistical errors associated to 
the measurement of the cumulants are negligible in the scales considered. 
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Figure 3. Ground-state energy per particle for the 2x3 FNU hard-core boson 
Hubbard model versus the interaction strength 7 with Np = 2 and A'p = 3 
particles per spin. We compare the results obtained by solving Eq. 1881 at 
truncation orders fcmax = 1, 2, 3, 4 (different lines) with those from exact numerical 
diagonalizations {-|-). The statistical errors associated to the cumulants, evaluated 
with A' = 200 and M = 10^, are negligible in this scale. 
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Figure 4. As in Fig. 13 for the 4x4 lattice with Np = 5 and Np = 8 particles 
per spin. We compare the results obtained by solving Eq. 1881 at truncation 
orders fcmax = 1,2,3,4 (different lines) with those from quantum Monte Carlo 
simulations ( X ). The statistical errors associated to the cumulants, evaluated with 
N = 200 and M = 10^ , and to the quantum Monte Carlo results are negligible in 
this scale. 



6. Conclusions 



By using saddle-point techniques and a cumulant expansion theorem, we have 
exploited an exact probabilistic representation of the quantum dynamics in a lattice 
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Figure 5. Enlargement of Fig. |21 (left panel) and Fig. |1] (right panel) in the 
small 7 region of the systems at half filling considered there. We show the results 
obtained by solving Eq. I88i at truncation orders femax = 2, 3, 4 (different lines, left 
panel) and fcmax = 2,3,4,5 (different lines, right panel) in comparison with those 
from exact numerical diagonalizations (+, left panel) and quantum Monte Carlo 
simulations (dots with error bars, right panel). The statistical errors associated 
to the cumulants are negligible. 



to evaluate the matrix elements of the evolution operator of a system of hard-core 
bosons in the Umit of long times. The approach yields a simple scalar equation for 
the ground-state energy in terms of the asymptotic cumulants S'-*^-' of the values of 
the potentials, V , and of the kinetic quantities, T, assumed by the system during 
its long-time evolution. Since the cumulants depend only on the structure of the 
system Hamiltonian, once they are known, this equation provides the ground-state 
energy and, via the Hellman-Feynman theorem lOl, the ground-state expectation of 
any other operator, analytically as a function of the Hamiltonian parameters. In 
contrast, quantum Monte Carlo methods require, due to the unavoidable branching 
or reconfiguration techniques (see [7j and and references therein), different 
simulations for different values of the parameters. 

The analytical character of the present approach suggests many potential 
applications. Here, we briefly envisage some of them. 

i) The knowledge of the ground-state energy as a function of the Hamiltonian 
parameters, Eqb = Eqb{^), allows in principle the determination of ground-state 
quenched-averages, / f{EoB{^))dP^, of crucial interest in the study of disordered 
systems. Here, is a given probability density of the disorder parameters ^. Note 
that when the dimension of the space of the Hamiltonian parameters is larger than 1, 
the evaluation of the above integral by a quantum Monte Carlo approach may be a 
hard, if not unmanageable, numerical task. 

ii) Even though we have specialized the study to the easiest case, namely that of 
hard-core bosons in the absence of a magnetic field, as evidenced at the beginning of 
Section HI and in Section HI A, our approach is general and not limited to hard-core 
boson systems. Fermions and bosons in a magnetic field may in principle be treated 
in a similar way, provided that we properly take into account also the multiplicities 
Nx, see Eq. (|18a|l . This possibility is of great interest since, as well known, fermions 
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and bosons in a magnetic field are both affected by the so called sign (or phase) 
problem, which in practice inhibits the accomplishment of unbiased quantum Monte 
Carlo simulations. 

Hi) We note that our approach is also analytical in the time parameter, t. By 
properly taking into account the derivatives of the cumulants with respect to the 
number of jumps, QatS^^'*^). one can obtain not only the asymptotic behavior of the 
matrix elements, Eq. 1)83(1 . but also their behavior at finite times, t, both imaginary 
or real. The latter possibility constitutes another chance of great interest because in 
general, due to the presence of oscillating terms, also the real time behavior is affected 
by a sort of sign problem and the quantum Monte Carlo simulations are reliable only 
for short times [7]. Of course, the complete knowledge of the time behavior would 
imply that of the excited states. 

iv) Finally, we mention a different possible application of the result ((88|l . This 
equation can be exploited in a Monte Carlo framework for effectively sampling the 
ground-state energy (work in progress). Essentially, the better efficiency of this Monte 
Carlo method, compared with that directly deduced from Eq. (|12|l . see Ref. [7] for 
details, follows from the fact that the stochastic times of the original probabilistic 
representation have been analytically integrated out, as done in Section III B, so that 
the fluctuations are necessarily reduced. 
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Appendix A. Proof of the cumulant summation rules 

In this Appendix we prove the relations ((60|l and H61|l . Due to the constraints H56|l . it 
is trivial to sec that the statistical moments (non-connected correlation functions) of 
order k, {vai ■ ■ ■ ^a^) cti, ■ • • , a/c G J^iC ^ satisfy the following summation rules 

for A: = 1, whereas for fc > 1 

Since (t'a)^^ = {va) Eq. (|A.ip coincides with Eq. 

To demonstrate Eq. ((61|1 . let us introduce the short notation 

m(/W) = (i.„,...z.„,)^ (A.3) 

s(/W) = (^„,...z.«ji^\ (A.4) 

where I^^'^ = {1, . . . , /c} is the set of integers that appear as subscripts in ai . . .Uk- 
More generally, for any nonempty subset Ip of Ip = {pi,p2, ■ ■ .}, we define m(Ip) 
and s(/p) as 

milp) = (^Va^^ Va^^ . . (A.5) 
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S{lp) = (v^p^Vap^ ...J^ . (A.6) 

The cumulants can be defined implicitly in terms of the moments according to the 
relation [5| 

Upip=iw p 

where the sum is extended to all the unordered decompositions of the set j'*^' in 
disjoint nonempty sets Ip such that Uplp = 

Let us proceed inductively and suppose that Eqs. (|61|l hold for any value of the 
order k such that 2 < k < n — 1, i.e. 

^s(/('=)) = 0, 2<fc<n-l, (A.8) 

Qfc 

where the sum runs either over the sets or 
At the order n, we rewrite Eq. I|A.7() as 

„(/(")) = .(/("))+ J2 U^iipMM) 

Up/p=/<"-l) P 

+ E n^(^p)' (A.9) 

where {n} is the set having only the element n and s({n}) = (j^c(„)^^ By summing 
over the index q;„, for a„ € or a„ S and using the relations ljA.2|) together with 

m(/("-i)) = ^. (/("))+ ^ l[sil,) 

Up/p=/("-l) p 

+ E nE^(^^')- (A.io) 

Up/p=/("), /p#/<"), Ip9^in} P 

According to Eq. (|A.7p . the second term in the r.h.s. of Eq. ljA.10|) is equal to 
m(/'"~^-') whereas the third term involves only sets Ip with < n — 1 so that 
by using the inductive hypothesis, Eqs. IjA.Sp . we have 

E^(^'"') = o- (A-ii) 

Finally, it is easy to verify by direct inspection that in the case k = 2 wc have 

^s(/(2)) = 0, (A.12) 

so that Eq. (|61|l is proved. 
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